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X-ray computed tomography (CT) has an indispensable role in constructing 3D images of ob¬ 
jects made from light materials. However, limited by absorption coefficients, X-rays cannot deeply 
penetrate materials such as copper and lead. Here we show via simulation that muon beams can 
provide high resolution tomographic images of dense objects and of structures within the interior of 
dense objects. The effects of resolution broadening from multiple scattering diminish with increas¬ 
ing muon momentum. As the momentum of the muon increases, the contrast of the image goes 
down and therefore requires higher resolution in the muon spectrometer to resolve the image. The 
variance of the measured muon momentum reaches a minimum and then increases with increasing 
muon momentum. The impact of the increase in variance is to require a higher integrated muon flux 
to reduce fluctuations. The flux requirements and level of contrast needed for high resolution muon 
computed tomography are well matched to the muons produced in the pion decay pipe at a neutrino 
beam facility and what can be achieved for momentum resolution in a muon spectrometer. Such an 
imaging system can be applied in archaeology, art history, engineering, material identification and 
whenever there is a need to image inside a transportable object constructed of dense materials. 


INTRODUCTION 


Computed tomographic (CT) images are formed by 
combining X-ray projection images from multiple angles, 
each of which record the amount of X-ray absorption 
as function of position, which is proportional to line- 
integrated X-ray absorption coefficient (J. However, the 
application of CT is limited by the large absorption 
coefficient of X-rays in certain materials such as heavy 
metals. For example, for 100 keV X-rays the attenuation 
coefficient is ~ 3.5 cm -1 [2]. After 2 cm of copper 
plate, the intensity is reduced to 0.1%. If the object 
is significantly thicker, then the required power of the 
X-ray generator increases exponentially [3]. 

The idea of using muons for imaging dense objects dates 
back to Luis Alvarez and collaborators, who used cosmic 
muons to search for hidden chambers in the ancient 
Egyptian pyramids 0]. Cosmic muons continue to be 
used for diverse applications from geological information 
from imaging volcanos [5, 6| to commercial and security 
use in the detection of dense radioactive materials 
in cargo containers [7H9]. Most cosmic muon imaging 
experiments are, however, counting experiments and 
the muon energy and direction are not well controlled. 
The cosmic muon flux is not suitable for imaging small 
objects with high resolution. Therefore, we would 
like to investigate the possibilities and limitations of 
using muon beams from an accelerator complex for the 
purpose of high resolution tomographic imaging. 

Muons incident on matter with electron number den¬ 
sity n, atomic number z, and mean ionization energy I 
will lose energy according to Bethe-Bloch equation [TO]: 
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Depending on the internal structure of the object, muons 
incident on different regions will traverse materials of dif¬ 
ferent types and densities, and as a result come out with 
different mean energies. If the mean energy is measured 
as a function of the incident position, we get a projected 
map of the muon energy loss. If this projection is taken 
from different directions, tomographic reconstruction al¬ 
gorithms, such as filtered back- project ion pQ, can be used 
to reconstruct the internal structure of the object. Since 
the reconstruction algorithms are well-studied and be¬ 
yond our scope, we will demonstrate the result without 
focusing on the reconstruction. 

It is worth pointing out that similar approach exists for 
protons [lTj, but nucleon-nucleon interaction will limit 
the range of protons for the energy and material of inter¬ 
est, whereas muons interact mainly electromagnetically 
and are free from such constraint. 


METHOD 


Following the conventions of CT, we use a modulation 
transfer function (MTF) [12] to characterize the perfor¬ 
mance of muon imaging. MTF measures how the details 
of different length scales are modulated. MTF is defined 
as the Fourier transform of the line-spread function pffi] . 
which describes how an infinitely sharp line will be dis¬ 
torted by the imaging system. The line-spread function 
is the derivative of the edge-spread function, which de¬ 
scribes how an infinitely sharp edge gets distorted. We 
will start with the edge-spread function. 

The test object is two rectangular layers of heavy metal 
(copper or lead) sandwiching a third layer, half of which 
is made of the same heavy metal, and the other half- 
filled with air (see FIG. [l]). On average muons incident 
on the heavy metal side lose more energy than the less 
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FIG. 1: The test object consists of an inner layer 
sandwiched between two outer layers. Half of the inner 
layer is made of the same material as the outer layer. 
Muon beams are scanned through the test object and 
the energy loss is measured by muon spectrometers. 


dense side, and an edge is formed. Normally a numeri¬ 
cal derivative should be taken to obtain the line-spread 
function[13]? but we found that the derivative will intro¬ 
duce a large error. So instead we fit the shape of the 
edge with an error function and the line-spread function 
is obtained by taking the derivative of the fitted error 
function. This approach is less prone to numerical er¬ 
rors while retaining the main features. In addition, the 
Fourier transform of the Gaussian line-spread function is 
a Gaussian MTF. In the remaining step, the standard de¬ 
viation of the Gaussian MTF is used as a figure of merit 
for the imaging resolution and will be referred to as the 
modulation transfer coefficient (MTC). 

In the Geant4 simulation [14], 6 x 10 6 muons are passed 
through the test object. Muons with the same source 
position (bin) are averaged and the average muon energy 
after transmission is plotted as function of source posi- 
tion (FIG.|3). 


RESULT 

An example edge-spread function and MTF are 
shown in FIG. [2] The corresponding projected images 
are shown in FIG. [3j From FIG. [2a] as the energy is 
increased, the edge becomes sharper. This is visible in 

fig. m 

The MTC as a function of the incoming muon energy 
is shown in FIG. [U The resolution increases with muon 
energy, and thin layers are better than thick layers. 
For the same thickness and muon energy, copper is 
better than lead. These findings are consistent with the 
expected effects of multiple scattering. The broadening 
is inversely proportional to the momentum of the muon 
and increases with the square root of the thickness of 
the material in units of radiation length flOl. 


Edge Spread Function 
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(a) Edge-spread functions. 
Modulation Transfer Function 



(b) Modulation transfer functions. 


FIG. 2: Edge-spread functions ( |2a| ) and MTFs (2bJ) for 
muon energies of 600 MeV, 1 GeV and 3 GeV. The test 
object material is copper with an air gap in the middle. 
The inner and outer layers are 10 cm thick. 


The relative amplitude, defined as the ratio of the 
amplitude of the fitted error function to the vertical 
offset, is shown as function of muon energy in FIG. [5) 
As the muon energy increases, the relative amplitude 
(contrast) decreases. The relative variance, defined as 
the variance with respect to the fitted error function 
divided by the amplitude increases after reaching a 
minimum. For high energy muons, image resolution is 
improved at the cost of requiring better detector resolu¬ 
tion to resolve the reduced contrast. Larger statistics are 
required to reduce fluctuations in a high resolution image. 

Sensitivity to density transition edges were tested for 
1 GeV muons, a 5 cm inner layer, 10 cm outer layers, 
and where the air gap is replaced with copper (lead) of 
different densities. The density is scaled with respect to 
the normal density of copper (lead). The MTC as a func¬ 
tion of density shown in FIG. [6] The MTCs are largely 
independent of the densities, except when the densities 
are the same to within 5%. The relative amplitude as a 
function of density fraction is shown in FIG. [ 7 ] When the 
imaging materials are of the same type but have differ¬ 
ent densities, the closer the densities are, the higher the 
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Average Muon Energy 



(a) 600 MeV muons. 

Average Muon Energy 



(b) 1 GeV muons. 

Average Muon Energy 



(c) 3 GeV muons. 

FIG. 3: Muon images with different energies incident on 
a test object with a 10 cm inner layer of air on the 
right-half and copper on the left-half and sandwiched 
between two outer layers of 10 cm copper plates. As the 
incident muon energy is increased, the boundary 
becomes sharper at the cost of reduced contrast. 


required detector resolution. For example, if we had a 
muon detector with a resolution of 2% at 1 GeV, we can 
resolve 8.9 g/cm 3 copper from copper scaled to a density 
of 7.12 g/cm 3 . 

In FIG. [8] we demonstrate a 3D reconstruction of a 
test object made of concentric spheres of lead, iron and 
copper, together with its cross-section image and radial 
distribution of muon absorption coefficient. Between dif¬ 
ferent layers there is 1 cm of water. Reconstruction was 
done with filtered back-projectionpQ with 200 angles. 


DISCUSSION 


The resolution is limited by the multiple scattering of 
muons inside the materials, so any measure that reduces 


Modulation Transfer Coefficient 



FIG. 4: The modulation transfer coefficients(MTC) are 
plotted as a function of muon energy for several test 
object configurations. The test object has a 5 cm inner 
layer made of copper (or lead) and air. The two outer 
layers are made of copper (or lead) and each layer is 
10 cm (or 20 cm) thick. The MTC increases linearly 
with incident muon energy and larger statistics are 
needed at higher energies. At the same energy, thin 
layers are better than thick layers, and for the same 
thickness copper is better than lead. 




Muon Incident Energy (GeV) 


FIG. 5: Relative amplitude (top) and relative variance 
(bottom) as function of incident muon energy. As the 
muon energy increases, the relative amplitude decreases 
while the relative variance first reaches a minimum and 
then increases. Low relative amplitude requires higher 
detector resolution, and large variance requires more 
statistics. 


muon scattering will improve the image resolution. Once 
the geometry is fixed, the only available option is to 
increase the muon energy [10]. However, this has a cost 
of reduced contrast and increased noise, and requires 
detectors with higher resolutions and a larger statistics 
of muons. This trade-off can benefit from a spectrum of 
muon momenta. One can scan the material with lower 
energy muons to determine the approximate regions 
and then switch to higher energy to determine the 
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FIG. 6: The top (bottom) curve shows MTC for 5 cm 
copper (lead) of different densities with 10 cm copper 
(lead) layers. The muons have energy 1 GeV. The 
density fraction refers to the scale to the normal density 
of respective materials. At a density fraction close to 1, 
the fluctuation grows since it becomes harder to 
distinguish the two different densities. While copper 
and lead have quite different MTCs, they stay constant 
for a large range of densities. 



FIG. 8: (left) 3D reconstruction of a test object using 
Visit [15], taken with 3 GeV muons. The object consists 
of 5 concentric spheres, with thickness 3, 9, 15, 21, and 
30 mm, respectively, and 10 mm of water between 
different metallic layers. Each metallic layer consists of 
lead, copper and iron of equal thickness, (right) 
Cross-section view of the test object, and the muon 
absorption along a radial line. For the outer layers, we 
can identify three distinct materials indicated by the 
color in the cross-section, and the steps in the line plot. 



FIG. 7: Contrast and variance as function of density 
fraction. As density fraction gets close to 1, the 
contrast drops and the variance increases. Since the 
MTC stays approximately constant, the ability to 
distinguish materials of the same type but different 
density is limited by the detector resolution. 


boundaries with better resolution. In these studies, no 
selection is made. Resolution can be further improved 
by discarding muons with large amount of scattering. 

Since muon tomography records the information of 
individual muons, one can use the amount of muon 
scattering to identify internal boundaries where variance 
come to local maximum, and to identify different ma¬ 
terials. The latter has been demonstrated with cosmic 
muons0 and can be applied to accelerator muons. In 
addition, for muons that are monochromatic or have a 


well-measured incident momentum, from the measured 
energy loss one can use he Bethe-Bloch equation (Eq. [l]) 
to identify materials or test models of different material 
compositions. In this way, one can examine the radiog¬ 
raphy of an object. 


Muon imaging systems can be co-located at neutrino 
beam facilities. In accelerator-based neutrino experi¬ 
ments, muon neutrinos are obtained from charged pion 
decays along with muons [16]. However, these muons 
are typically filtered from the beam without being used. 
Given the accelerator neutrino energy, the muon energy 
can be calculated. Let 0 denote the angle between neu¬ 
trino and the pion and assuming zero neutrino mass, neu¬ 
trino energy is: 

o 9 

E = 

2E n — 2p n cos 6 
(1 — mj'/rni) 

~ — 1 - ^ 02 — E *- ( sma11 angle) 

A similar calculation is done in m- In the center-of- 
mass frame, the neutrino angular distribution is isotropic, 
so in the lab frame the neutrino direction has angular 
distribution P(6) = sin0/[y 2 (l — f3cos6) 2 ], with the most 
likely angle for neutrino 0 * = cos -1 /3. The ratio of muon 
energy to the neutrino energy for the most probable angle 
is: 


( 2 ) 

( 3 ) 


E/J,/Ey 


2 < 

m\ — m 


2 

A 6 


— 1 


3.7. 


( 4 ) 


Taking NuMI at Fermilab as an example, the neutrino 
beam energy is in the range 1 — 3 GeV for the low 
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energy option[I8j [19]. The corresponding muon energy 
is 4 — 10 GeV, which lies in the energy range of a muon 
computed tomography system. The CNGS at CERN has 
about 3 times larger energy 20.121]. The required muon 
energy is different for each object and material, but 10 
GeV is enough for imaging most objects of moderate 
size. If there is a need to reduce the muon energies 
coming out of the accelerator facility, then the muon 
beam can be cooled down with an absorber and selected 
to the desired energy with a spectrometer. 

The muon flux produced in associated with the pro¬ 
duction of neutrino beams at the NuMI facility is roughly 
10 7 /cm 2 per spill, over a roughly 1 m 2 beam area pTHl8] . 
The NuMI spill cycle is a spill duration of approxi¬ 
mately 10 /is every 1.87 s[T8], In between each spill, 
the orientation of the object with respect to the beam 
can be stepped in angular increments. For one degree 
steps in azimuthal and polar directions, the imaging data 
for a roughly 1 m 2 cross-sectional area can be collected 
in roughly 12 minutes. The test object reconstruction 
in FIG. [8] assumes approximately 20 min beam time. 
The limitation of such a device is largely instrumen¬ 
tal. To simultaneously acquire such a large flux of muon 
data requires highly pixelated silicon tracking operat¬ 
ing at hundreds of MHz. There are silicon pixel read¬ 
out systems with these capabilities designed for high-rate 
environments[22]. Spectrometer selection could aid in re¬ 
ducing the flux by selecting a narrow band of incident 
muon momenta. 


CONCLUSION 

In this paper, we have shown that computed tomog¬ 
raphy with accelerator muons can be used in place of 
X-rays to overcome the limitations of imaging inside 
of objects made of heavy metals such as copper and 
lead. The spatial resolution is shown to increase with 
increasing muon energy, due to the reduction in multiple 
scattering. However, this reduction is at the cost of 
image contrast and requires more statistics to suppress 
fluctuations. The resolution also depends on the geom¬ 
etry, material and material boundaries. In the case of 
copper and lead, lead scatters muon more, and therefore 
results in poorer resolution. For materials of the same 
type but different densities, the resolution remains some¬ 
what constant as a function of the density fraction. The 
relative amplitude decreases as the densities get closer, 
and, therefore, places a more stringent requirement on 
the detector resolution to distinguish regions having 
different densities. Either by measuring the multiple 
scattering or energy loss of individual muons, it is 
possible to do material identification. 

We have framed the capabilities of muon computed 


tomography through simulation. As with X-ray CT at 
the outset, it takes foresight and planning to envision 
the breakthroughs that dense object imaging may yield. 
Muon CT is a probe like no other into the unknown that 
hides deep within dense structures. The muons produced 
in association with high intensity neutrino beams fall into 
the energy range of interest for muon imaging. This pro¬ 
vides a unique opportunity to incorporate muon imaging 
systems into existing or future neutrino beam facilities. 
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